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Abstract 

We compute the full non-perturbative ghost and gluon two-point Green functions by using gauge 
field configurations with Nf = 2 and Nf = 2+1 + 1 twisted-mass quark flavours. We use simulations 
with several different light quark masses, heavy quark masses close to that of the strange and 
charm quarks, and lightest pseudoscalar masses ranging from 270 to 510 [MeV]. Quark flavour 
effects on both the gluon and the ghost propagators are then investigated in a wide range of 
momenta, bridging the deep infrared and intermediate momenta domain of QCD interactions in 
the presence of dynamical quarks. The ghost-gluon vertex is also indirectly probed through a 
consistency requirement among the lattice data for the gluon and ghost propagators and the ghost 
propagator Schwinger-Dyson equation. The effective full QCD coupling is finally constructed, and 
its dependence on the presence of dynamical fermions scrutinized. 



PACS numbers: 11.15.Tk, 12.38. Gc, 11.15.Ha, 12.38.Aw, 14.70.Dj 



I. INTRODUCTION 



During the past few years, lattice simulations have considerably improved our understand- 
ing of the infrared (IR) sector of non-abelian Yang Mills theories. In particular, quenched 
Landau gauge simulations [IH5], performed on lattices with large volumes, have unequivo- 
cally demonstrated that the gluon propagator saturates in the (deep) IR region. This is true 
for the space-time dimensions d = 3 as well as 4, irrespectively of the number of colors N c of 
the gauge group SU (N c ) under consideration. At the same time, the ghost dressing function 
effectively acquires its tree-level behaviour, with the functional form of the propagator being 
~ l/g 2 . 

Within the continuum formulation of the theory, these lattice results are in agreement 
with the solutions of the corresponding all-order Schwinger-Dyson equations (SDEs) [6], [7j 
and exact renormalization group (RG) equations [8]. Other approaches such as the so- 
called refined Gribov-Zwanziger formalism j9] also converge to the same conclusions. This 
has caused a paradigmatic shift among practitioners: the gluon is now thought to acquire 
a momentum-dependent mass m(q 2 ) whose magnitude can be large at IR momenta, but 
vanishes with increasing spacelike momenta (i.e., q 2 ^> Aq CD ), thereby maintaining full 
accord with perturbative QCD. Gluon confinement is then realized, as it is customarily 
done in the case of quarks, through the violation of reflection positivity (signaled by the 
presence of an inflection point of the propagator scalar cofactor A(g 2 )) instead of achieving 
an area law for a Wilson loop or a linearly rising potential (criteria which are irrelevant 
to the question of light-quark confinement [H)]), or satisfying ad-hoc criteria involving the 
ghost sector (which, as already pointed out above, completely decouples in this regime). 

The extension of these quenched lattice results to full QCD, i.e., to a non-Abelian SU(3) 
theory with the inclusion of dynamical quarks, has not been extensively pursued, neither 
in the continuum nor on the lattice. In the former case, a first analysis of the effects on 
the gluon propagator due to dynamical quarks has recently been reported in [TT] within 
the so-called PT-BFM (pinch technique-background field method) truncation scheme |12l - 
[15]. Earlier related endeavours on the lattice can be traced back to [16] where an 0(a 2 ) 
Symanzik-improved action with 2 + 1 staggered fermion flavours was employed, and [TT] 
where a tadpole-improved gauge action with 2 dynamical overlap fermions was used instead. 
However, an independent affirmation of these results by implementing different lattice ac- 
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tions 1 as well as their extension for different numbers of flavours, has been a pending issue 
since then. 

This article provides a comprehensive quantitative study of the aforementioned Green 
functions which incorporate the effects stemming from the presence of dynamical quarks. 
To this end, we compute the gluon and ghost two-point Green functions from the gauge 
configurations generated by the ETM collaboration [T9l 120] for the cases of (i) two light 
degenerate quarks (Nf = 2) and (ii) two light and two heavy 2 (Nf = 2 + 1 + 1) mass- 
twisted lattice flavours [21] • Furthermore, we apply our lattice results to carry out an 
indirect study of the ghost-gluon form factor (as done for quenched lattice data in [22]), 
by employing a hybrid approach where the solutions of the ghost SDE are studied using 
the gluon propagator determined in our simulations as an input. Consequently, the natural 
requirement to reproduce the lattice ghost dressing function data from the corresponding 
SDE solution will pin down the ghost-gluon vertex form factor, which will be shown to 
deviate considerably from its tree-level value. The constructed SDE solutions then allow us 
to extrapolate the lattice ghost data down to the vanishing momentum region and obtain 
reliable information on the saturation point of both the ghost dressing function as well as of 
the so-called Kugo-Ojima parameter [23]. Finally, the QCD effective charge, defined in [24J, 
is computed by properly combining the gluon propagator and the ghost dressing function 
with the lattice estimate of the coupling in the so-called Taylor scheme (e.g., see [25]) at a 
given (large enough) momentum. 

The main results of this article can be summarized as follows: 

• The effect of the presence of dynamical quarks on the gluon propagator A is twofold: a 
suppression of both the "swelling" region at intermediate momenta and the saturation 
value in the deep IR (which can be interpreted as the gluon becoming more massive in 
the presence of quarks). In addition, one observes that the more light flavours there 
are, the bigger the effect is, which is in accordance with what we would naturally 
expect. Light virtual quarks can be copiously produced, thus screening the interaction 

1 Some preliminary results obtained from simulations with large lattice sizes (far from the continuum limit) 
and Nf = 2 Wilson-Clover fermions have also been reported in [15] . 

2 It should be also noticed that these 2 + 1 + 1 configurations provide a realistic simulation of QCD below 
the bottom quark mass threshold, mainly at the momentum scales which we compute the Green functions 
for. 
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and suppressing the very same mechanism which triggers gluon mass generation. As 
the fermion mass is increased (at a fixed flavour number) the effect gets smaller, since 
the heavier the fermions, the lesser is the statistical likelihood of their pair-production. 
At a sufficiently large value of their mass, they essentially decouple and the gluon mass 
generation is practically insensitive to their presence. With respect to this point it 
should be noticed that our results turn out to be in agreement with the SDE study 
reported in [11] confirming at the same time the general trend reported in the earlier 
lattice studies of 



• On the other hand, the effect on the ghost dressing function F is much milder and is 
diametrically opposed to the one encountered for the gluon case, i.e., it consists in a 
small increase of the saturation point. This result is also in harmony with what one 
would intuitively anticipate. In the SDE for the ghost, the quark propagator does not 
enter directly, but only It does so only through the gluon propagator or via higher loop 
corrections to the gluon-ghost vertex. Therefore, it is natural to expect the influence 
of dynamical quarks to be less pronounced for the ghosts. 

• When the gluon propagator obtained is used as an input in the ghost SDE, one finds 
that the requirement for the SDE solution to match the ghost propagator lattice data 
naturally provides a stringent check on the ghost-gluon vertex; specifically, this exercise 
will show that this vertex differs significantly from its tree-level value. 

• Finally, when all the results are used to form the RG invariant combination aAF 2 
eventually leading to the QCD effective charge, we observe that, although obviously 
modifying the ultraviolet (UV) parameters controlling the running of the coupling and 
its magnitude, the number of fermions flavours does not affect the IR behaviour of 
this quantity. 

The paper is organized as follows: Section [TTJprovides the reader with some of the technical 
details of the lattice set-up used for the computation of the relevant gluon and ghost Green 



functions. Next, in Section III, we present the results of the simulations, emphasizing the 



differences with respect to the quenched results; volume artifacts are also addressed in some 



detail. The ghost SDE is then solved in Section IV, and the effective coupling evaluated in 



Section [V] Finally, we provide the conclusions in Section VI 
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II. GENERALITIES 



The following section is a reminder of how the ghost and gluon propagators are computed 
from the lattice simulations of gauge fields for light and heavy mass-twisted lattice flavours. 
It should be noticed that these propagators have been obtained (but not presented) earlier, 
as a by-product of the computation of the running coupling in the momentum subtraction 
(MOM) Taylor scheme [26H28] . These references, which the interested reader is referred 
to, also contain relevant details concerning lattice actions, set-ups and the treatment of 
artifacts. 

In our simulations, the lattice fermion action for the doublet of light degenerate quarks 
is given by [29J 

Si = a 4 ^2xi(%) (D W [U) + m j + ina 5 T 3 )xi(x), (2.1) 

X 

whereas, for the heavy doublet, we employ 

S h = a 4 ^2x h (x) (D W [U) + m 0)h + i^alsn + /i 5 r 3 ) Xh{x), (2.2) 

X 

where Dpj/[C7] stands for the standard massless Wilson Dirac operator. In the gauge sector, 
the tree-level Symanzik improved gauge action (tlSym) [30J is applied for Nf = 2 and the 



Iwasaki improved action [3T1 132] for Nf = 2 + 1 + 1. In addition to the plaquette term U^ u . 



rlxl 

X,/J,,U 

this formulation of the action also requires including rectangular (1x2) Wilson loops U^ 2 U . 
For instance, in the tlSym case, the action reads 

s » = f E \ b ° E i 1 - ReTr (tOl + 6 i E t 1 - ReTr (Cy] }> (2-3) 

l<H<v fij^u 

where = g is the bare lattice coupling and one sets b\ = —1/12 and b = 1 — 8&i as 
dictated by the requirement of continuum limit normalization. Configurations of the gauge 
fields generated by the above actions are next gauge fixed to the (minimal) Landau gauge. 
This is done through the minimization of the following functional [of the SU(3) matrices 



Re 



1 ~ ^{^U^g^x + fi) 



(2.4) 



with respect to the gauge group element g. 

To get as close as possible to the global minimum, we apply a combination of an over- 
relaxation algorithm and Fourier acceleration, considering the gauge to be fixed when the 



condition |<9 M A M | 2 < 1CT 11 is fulfilled and the spatial integral of Aq is constant in time to 
better than 10~ 6 . Evidently, this procedure cannot avoid the possibility that lattice Gribov 
copies are present in the ensemble of gauge fixed configurations. However, extensive litera- 
ture in the quenched case (see for example pi]) shows that such copies do not seriously affect 
the qualitative and quantitative behavior of the Green functions in question. Given also the 
relative large physical volumes simulated, we will proceed under the working assumption 
that this feature survives unquenching, as was also verified in [18]. 

After the lattice configurations have been projected onto the Landau gauge, one can start 
calculating the Green functions of interest. 

To begin with, we consider the gluon propagator. The gauge field is defined as 

UJx) - Utix) 1 UJx) - UJx) 
A,(x + m = - ' " - tt Tr ^ J 2.5 

2iag 3 2tag 

with fx indicating the unit lattice vector in the \x direction. The two-point gluon Green 
function is then computed in momentum space through the following Monte-Carlo average 

Aj£(ff) = (A«(q)Al(-q)) = 6 ab (V - q -^\ A(g 2 ), (2.6) 

with 

A«(q) = l - Tr M* + A/2) exp[zg ■ (x + fi/2)]X a . (2.7) 

X 

In the formula above A a are the Gell-Mann matrices and the trace is evaluated in color space. 

The Landau gauge ghost propagator can also be computed in terms of Monte- Carlo 
averages of the inverse of the Faddeev-Popov operator, i.e., 



i <Vexp[- <- ..*f*,-i\A-x* F tf 
with M written lattice divergence 



^V) = ^ ( £ expfa • (x - y)] (M~X y ) = ^^f 1 , (2-8) 

\ x,y I 



M{U) = -^V-D{U) ) (2.9) 

and the operator D acting on an arbitrary element of the Lie algebra rj according to 

D(U)rj(x) = l - [U.ixMx + fj,)- ri(x)U„(x) + V (x + ^ - Ul(x) V (x)] . (2.10) 

More details on the lattice procedure for the inversion of the Faddeev-Popov operator can 
be found in 1 33 1 - 
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Next, if we indicate with A the regularization cutoff (e.g., A = a _1 (/3) if one specializes to 
lattice regularization), one can obtain the renormalized gluon propagator and ghost dressing 
function as 

A R (gV) = lim Z 3 - 1 (/ i 2 ,A 2 )A( g 2 ,A 2 ), 

A— >OD 

F K (q\^) = lim Z 3 -V,A 2 )F(g 2 ,A 2 ), (2.11) 

A— K» 

where one imposes the standard MOM renormalization conditions 

Ah.Gu 2 ,/! 2 ) = 1/ M 2 ; F R (/iV) = 1. (2.12) 

When unnecessary, we will refrain from explicitly indicating the renormalization point de- 
pendence of the various renormalized quantities. 

We conclude this section by commenting briefly on the crucial role played by the so- 
called if(4)-extrapolation procedure |MH3S], which have been used to correct the data for 
discretization artifacts (otherwise plaguing the reliable determination of A and F) due to the 
breaking of the 0(4) rotational invariance down to the H(4) isometry group. Specifically, let 
us observe that the gluon and ghost dressing functions (q 2 A and F) are dimensionless cor- 
relation functions, and therefore general dimensional analysis shows that they must depend 
on the (dimensionless) lattice momentum a q^, where 

271" 77 

Qfj, — ~~aj~^ > 7^ = 0,1,...,^, (2.13) 

being the number of lattice sites in the /x direction (in our case, N x = N y = N z = N t /2). 
However, if one considers a dimensionless correlator Q evaluated on the lattice, since 0(4) 
is broken down to H(4), one has 



g.- (a Y,a 2 ^,---) = g latt (aV)+ ° Q 



Q 



a 1 ' 

a 2 V+"--> (2-14) 



where = q^ is the first if(4)-invariant (and the only one relevant in the ensuing 
analysis). The if (4)-extrapolation procedure is thought to account properly for the breaking 
of 0(4) down to H(A) and thus recover the continuum-limit 0(4)-invariant result by means 
of the following prescription: one first averages over any combination of momenta being 
invariant under if (4) (a so-called H(4) orbit); next, one extrapolates the results towards 



the continuum limit (where the effect of a 2 q^ must vanish) by applying Eq. (2.14) to all 
the orbits sharing the same value of q 2 . The only assumption employed is that the slope 



coefficient in Eq. (2.14) depends smoothly on a q . 



/3 


^crit 


am 




CLflS 


{L/af x T/a 


confs. 


3.90 


0.161856 


0.004 






24 3 x 64 


50 


4.20 


0.154073 


0.002 






48 3 x 96 


50 


1.95 


0.161240 


0.0035 


0.135 


0.170 


48 3 x 96 


40 


1.90 


0.163270 


0.0040 


0.150 


0.190 


32 3 x 64 


50 



TABLE I: Lattice set-up parameters for the ensembles we used in this paper: k ci h is the critical 
value for the standard hopping parameter for the bare untwisted mass; ^ stands for the twisted 
mass for the two degenerated light quarks, while [i a and f^s define the heavy quarks twisted masses; 
the last column indicates the number of gauge field configurations we used. 



III. SIMULATION RESULTS 



In this section we describe the outcome of our lattice simulations. The parameters used 
are reported in Table [TJ The physical scale, i.e., the lattice size at any bare coupling /3, has 
been fixed by European Twisted Mass Collaboration (ETMC) through chiral fits to lattice 
pseudoscalar masses and decay constants. At the physical point, these are required to take 
on the values of f n and m n provided by experiments. The bare untwisted mass is tuned to its 
critical value by setting the so-called untwisted Partially Conserved Axial Current (PCAC) 
mass to zero, so that the twisted-mass fermions are at maximal twist. The renormalized 
running masses for light and heavy quarks are obtained from the bare twisted-mass as 

am 



V>c/s{<1q) = I7m~^ 7T\ ( a ^ =•= ry (_ y ^S ) > (3-1) 



a(/3)Z P (q ) ' 

a(p)Zp(q ) \ a Z P (q ) 



where go is the renormalization scale. The determination of the nonperturbative renormal- 
ization constants, in particular Zp and Zg, is the subject of an exhaustive computation 
program within the framework of ETMC (see for instance [37J for the Nf = 2 case and 
[55| I3~5] which contain some preliminary results for the Nf = 2 + 1 + 1 case). The degenerate 
light quark masses we used for the simulations (Table [T]), range from 20 to 50 [MeV], while 
the strange quark is roughly set to 95 [MeV] and the heavy charm to 1.51 [GeV] (in MS 
at go — 2 [GeV]). The lightest pseudoscalar masses for the simulations of Table [T] range 
approximately from 270 to 510 [MeV]. The biggest volume simulated corresponds to an 
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asymmetrical box of roughly 3 3 x 6 [fm 4 ]. 
A. Gluon sector 

The results obtained for the gluon propagator and dressing function for the cases of two 
light quarks and two light plus two heavy quarks are plotted 3 in Fig. [I] As far as the 
gluon propagator is concerned (top panel) one can clearly see the IR flattening typical of the 
massive solutions. However, when compared to the quenched case (shown for reference by 
the diamond-shaped gray data points), the propagator shows a less pronounced "swelling" 
at intermediate momenta and a lower freezing out value. 

To check the dependence of this latter effect on the lattice volume, we plot in Fig. [2] the 
value of A R (0) as a function of the inverse of the volume. Though we do not have enough 
simulations on large volume lattices to attempt any continuum extrapolation, it is evident 
that residual volume effects are expected to be small when the appropriate simulations (i.e., 
(3 = 4.20 for N f = 2 and both (3 = 1.95 and (3 = 1.90 for N f = 2 + 1 + 1) are considered. 
Furthermore, apart from the zero-momentum gluon propagator, the results for our two 
simulations in both cases appear clearly superimposed in the plots of Fig. [TJ indicating that 
volume effects are indeed under control. 

In addition, the quenched simulation can be viewed as an unquenched counterpart in the 
limit of infinitely massive fermions, and the Nf = 2 results as the limit of the Nf = 2 + 1 + 1 
case in the infinite mass limit of the heavy sector. Thus, one can unambiguously conclude 
that the presence of dynamical fermions suppresses the IR saturation point, and renders the 
gluon heavier. Also notice that the suppression tends to subside as the dynamical fermion 
mass increases. The decoupling of heavy fermions has been explicitly shown in the continuum 
through the SDE analysis of [TT], where it was found that the gluon propagator results for 
Nf = 2 + 1 approach those for Nf = 2 as the mass of the heavy flavour is increased (see Fig. 
17 of HU). 

Finally, the concave shape of the propagator ensures the violation of reflection positivity, 
thus implying that the unquenched gluon is also a confined excitation. 

The behavior of the dressing function (bottom panel) is similar. In this case, the greater 

3 If not stated otherwise we will be setting the renormalization point to be /i = 4.3 GeV 
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FIG. 1: The unquenched gluon propagator (top panel) and dressing function (bottom panel) for 
Nj = 2 (two light quarks) and Nf = 2 + 1 + 1 (two light and two heavy quarks). 

the number of dynamical quarks, the less pronounced the peak at the intermediate momenta. 
Analogously, the heavier the quark, the less the effect it entails on the overall shape of the 
dressing function. These results are in agreement with the SDE study of [11], as well as 
with the lattice findings of [T6HT8]. 
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H=4.3 [GeV] N,=2, p=4.05, V=24,x48 
N f =2, 0=3.90, V=24^x48 
N f =2, 0=4.20, V=48,x96 
N,=2+1+1, 0=1.90, V=32,x64 
N,=2+1 +1, 0=1 .95, V=48 J x96 



0.01 

l/V [fm- 4 ] 



0.1 



FIG. 2: The volume dependence of the IR saturation point of the gluon propagator A H (0) in our 
simulations. For the Nj = 2 case, we include an extra point corresponding to a simulation on a 
24 3 x 48 lattice, at /3 = 4.05 (k = 0.157010 and am = 0.006). 

B. Ghost sector 



The results for the ghost dressing function are plotted in the top panel of Fig. [3] In 
analogy with the quenched case, the data do not support a power-like singular behaviour in 
the (deep) IR region; rather one finds the characteristic freezing out feature of the massive 
solutions [3 HQ] . As one would expect on the basis of a naive perturbative analysis (there is 
no tree level coupling between ghosts and fermions) the effect of dynamical quarks on the 
ghost sector is much milder as compared to the gluon sector. 

The ghost dressing function F can provide valuable information with respect to the so- 
called Kugo-Ojima function [23]. This is due to a powerful identity dictated by the under- 
lying Becchi-Rouet-Stora-Tyutin (BRST) symmetry present in the continuum formulation 
of the theory, which leads to the relation (HI H2] 

F-\q 2 ) = l + G(q 2 ) + L(q 2 ), (3.2) 

where G(q 2 ) and L(q 2 ) are the form factors of a particular Green function A^ u (q) that plays 
a special role in the aformentioned PT-BFM truncation scheme jS], with 



KM = 6, u G{q 2 ) + q -^L(q 2 ). 



(3.3) 
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0.01 0.1 1 10 



q 2 [GeV 2 ] 

FIG. 3: The unquenched ghost dressing function (top panel) and the (approximate) Kugo-Ojima 
function (bottom panel) for Nf = 2 (two light quarks) and Nf = 2 + 1 + 1 (two light and two heavy 
quarks) . 

The important point here is that G(q 2 ) coincides (in the Landau gauge) with the Kugo- 
Ojima function |?TJ H2]. In addition, a detailed analysis of the L(q 2 ) form factor in the 
quenched approximation reveals that it is numerically subdominant in the whole range of 
momenta when compared to G(q 2 ) [12], and, furthermore, L(0) = 0. Since quark effects 
on A MJy (g) are suppressed, either due to their indirect presence in the full gluon and ghost 
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propagators, or in higher order corrections to the ghost-gluon kernel (the first one happening 
at the three- loop level in the kernel skeleton expansion, and therefore at four loops in A MI/ ), 
one naturally expects the same results to survive in the unquenched case, thus leaving us 
with the approximate relation 

G{q 2 ) ^ F~\q 2 ) - 1. (3.4) 

In the bottom panel of Fig. [3] we plot the function —G(q 2 ) and observe that its value at 
origin is practically unchanged when varying the number of flavours. Clearly the behavior 
is not dissimilar from the one revealed in quenched simulation, and the (extrapolated) IR 
saturation value looks once again far from the critical value 1 predicted by the scaling type 
solutions of the SDE and the related Kugo-Ojima confinement criterion. We will return to 
this issue in the next section. 



IV. GHOST SDE ANALYSIS 



In this section we carry out a hybrid analysis combining our lattice simulation results 
with SDE techniques, in a spirit analogous to what has been reported in [10] • The aim is 
to study the ghost sector in greater detail and, in particular, gain access to the ghost-gluon 
vertex form factor(s). As a welcome byproduct, we will obtain a reliable extrapolation of 
the ghost lattice data to the deep IR. 

Specifically, let us start by considering the ghost SDE, which can be recast in the following 
bare form 

1 d A k F(k 2 )A((k-q) 2 ) \(q-k) 2 



F{q 2 



(2*) 



k 2 (k — g) 2 



- k 



Q 1 



Hi(k,q), 



(4.1) 



where Hi(k,q) is non-longitudinal form factor of the ghost-gluon vertex, parameterized as 

ff c (-k,q;k-q) = ig Q f ahc k v ,f v , v {-k,q-k-q) 

= ig f abc [k v H 1 {k,q) + (k-q)„H 2 (k,q)), (4.2) 

with k and q being the outgoing and incoming ghost momenta respectively, and go the 
bare coupling constant. As explained in depth in [301 HH US], one can first renormalize the 



ghost and gluon propagators in Eq. (4.1), by using Eq. (2.12), and then apply a subtraction 



procedure to deal with the UV singularity of the ghost self-energy integral to obtain 
1 



F R (q 2 ) 



1 + z 2 z, g2 ° 



3^3— l t i dkK(k,q)Hr c (k,q)F R (k' 



(4.3) 
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where 



K(k,q) 



sin 4 #d# 



A R ((k-q) 2 ) A R ((k-pf 
{k-qf (k-p) 2 



(4.4) 



The renormalization point, fi 2 , is implicitly present as an argument for all the renormal- 



ized quantities. In obtaining Eq. (4.3), the subtraction procedure is applied for Eq. (4.1) 
evaluated at the two momenta k and p, both being parallel and such that p 2 = fi 2 . Hi in 



Eq. (4.3) is a bare but finite jl6] quantity which needs no renormalization while, in front 
of the integral, the renormalization constants and the bare coupling especially appear in 
the right combination to cancel the cut-off dependence and give the MOM Taylor scheme 
coupling (see e.g., [2"B"]). 



a T (/x 2 



9o - V 2, ,,2, 7 ,,2\ 
47T 



(4.5) 



This coupling a T (/i 2 ) for Nf = 0, 2 and 2 + 1 + 1 can be determined from lattice data 
(see, e.g., [25T - T28] ). In order to solve the ghost SDE in isolation (i.e., without coupling it 
to the much more complicated gluon SDE), one can use the just determined lattice gluon 



propagator A R as an input for the equation, thus fully determining the kernel (4.4). Now the 



only unknown term present in the equation is the ghost form factor Hi ; clearly the solutions 
to the ghost SDE will describe the lattice data with a better or worse agreement depending 
on our ability to model this form factor [221 E] ■ 

Through the analysis of the quenched lattice data, it was shown in [3Q] that the solutions 



of Eq. (4.3) grossly underestimate (by a factor of at least 2) the lattice data if one uses 



the tree-level value Hi = 1 for the ghost-gluon form factor. A constant does indeed do a 
better job [ID] but does not allow for a precise description of the deep IR behavior of the 
function [22] . This implies that a good description of the (quenched) ghost dressing lattice 
data calls for a ghost-gluon form factor with a non-trivial kinematical structure. Using the 
knowledge derived from the OPE analysis of |22j, coupled with the current lattice data on 
the (Landau gauge) ghost-gluon vertex [48], one can parameterize this form factor as 4 



Hi{k,0) = H° 



1 + 



N c g 2 {A 2 } k 2 
4(iV2 - 1) fc 4 + m 



(1 - Hi) 



w 



w A + k 4 



(4.6) 



Estimates for the gluon condensate 5 g 2 {A 2 ), and the IR mass scale miR, can be obtained 



Eq. (4.3) requires to be solved with the full vertex Hi(k,q), which is modelled in [22]; however it can be 



shown that H\(k, q) ~ H\(k, 0) is a good approximation to obtain the ghost dressing in the IR momentum 
region |4"5] . 
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H=3.61 [GeV] Kl f =0, p=5.70, various ^s 
N,=2, p=3.90, V=24,x48 
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N,=2+1 +1 ,_p=1 .95, V=48 3 x96 
SDE, full vertex (N,=0) 
SDE, full vertex (N,=2) 
SDE, full vertex (N f =4) 
SDE, trivial vertex (N f =0) 
SDE, trivial vertex (N f =2) 
SDE, trivial vertex (N f =4) 




0.1 1 
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FIG. 4: (Top panel) The ghost dressing function obtained from the solution of the ghost SDE (4.3) 
with the lattice gluon propagator as an input and a T = 0.25, 0.32, 0.37, respectively, for Nf = 0, 2, 
and 2 + 1 + 1 at \x = 3.61 [GeV]. Dashed lines correspond to solutions for the tree-level Hi(q, 0) = 1, 



while continuous lines to the inclusion of the full form factor (4.6). The latter is also shown in the 
bottom panel. The values of the parameters used to integrate the ghost SDE are : g 2 (A 2 } = 7 
GeV 2 , miR = 1.3 GeV and w = 0.65 GeV (the same IR ones for the three cases). Moreover, 
Hi = 1.26, 1.20, 1.18 for N f = 0, 2 and 2 + 1 + 1, respectively. 
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FIG. 5: Extrapolation of the (approximate) Kugo-Ojima function in the deep IR. The ghost 
dressing function F employed in this plot is generated by solving the ghost SDE. 

from lattice data and OPE analysis [17]; the constants H® and w (introduced in order to 
guarantee that Hi(0, 0) = 1, as suggested by current lattice data [IS]) can be adjusted so 



that the solutions Eq. (4.3) match the corresponding lattice data as closely as possible. 



The solutions of the ghost SDE (4.3) following the procedure just illustrated are pre- 
sented in Fig. [4] In the top panel of the figure one can clearly see that, similarly to the 
quenched case, a tree-level value for Hi does not give solutions which can describe the lattice 



data. However, once the kinematically non-trivial expression, Eq. (4.6) (bottom panel of 
the same figure), is included in the equation, one obtains an excellent agreement with the 
data. Therefore, effectively, the curves for Hi(q,0) represent a genuine prediction of our 
analysis. It would be interesting to confirm or refute this prediction through direct lattice 
calculations of the ghost-gluon three-point function. 

The good agreement between the SDE solutions and the lattice data allows for an extrap- 
olation of the latter towards the deep IR region, where one observes a very small increment 
of the saturation point (monotonic with Nf). This is particularly useful when scrutinizing 



5 The very notion of condensate have been recently questioned in [SOI [SI]- In particular, according to 
the new perspective suggested there, our gluon condensate g 2 {A 2 ) should be understood as a mass-scale 
parameter related to the local operator A 2 in the OPE expansion of the gluon Green functions. 
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the Kugo-Ojima function as shown in Fig. [5j which clearly depicts that the saturation point 
of the function is practically insensitive to the inclusion of dynamical fermions. 



V. EFFECTIVE COUPLING 

The results obtained for the gluon and ghost two-point functions allow us to extract 
the running of the full QCD effective charge for a wide range of physical momenta, and in 
particular in the deep IR region which is evidently inaccessible to perturbation theory. 

To begin with, let us recall that the QCD effective charge is defined, among practitioners, 
in primarily two different ways: the first one (to be denoted by a PT ) is obtained within the 
framework of the pinch technique [T5l 152] and represents the most direct generalization of 
the familiar QED effective charge concept to a non-Abelian setting; the second one (to be 
denoted by a gh ) corresponds to the non-perturbative generalization of the strong coupling 
in the Taylor scheme mentioned before [26T - T28] . 

The construction of either effective charges proceeds through the identification of a suit- 
able RG invariant combination. Before identifying this quantity however, let us observe 
that though the effective couplings a PT and a g h have a rather distinct theoretical origin and 
status, it turns out that, in the Landau gauge, they are related through the equation 



a gh (q 2 



1 + 



L(q 



2\ 



-2 



« PT (g 2 ); (5.1) 



1 + G(g*) 

evidently, in the approximation L(q 2 ) « 0, used throughout this paper, the two definitions 
coincide, and one has a PT (q 2 ) = a^q 2 ) = a(q 2 ). This implies also that one can choose as 
the RG invariant combination 

r(q 2 ) = a T (^)A R (q 2 ,^)Fl(q 2 ,^), (5.2) 



which can be readily obtained from the data presented so far [a T is given in Eq. (4.5)]. 
The quantity r(q 2 ) defined above is constructed in Fig. [6] for different number of flavours 
Nf] notice that for calculating the freezing out point r(0) the value of -Fr(O) has been 
extrapolated from the SDE results for the ghost dressing obtained in the previous section. 

A most salient feature of this plot is the absence (within the errors) of any flavour depen- 
dence in the IR region, more precisely starting from q 2 < 1 GeV 2 . Indeed, one observes that 
the flavour effects which control the behavior of the UV parameters of the theory (e.g., the 
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FIG. 6: The RG invariant quantity r(q 2 ) defined in Eq. (5.2) for the different values of flavours 
Nf. Errors for the Nf = case are underestimated, since we have used the ghost dressing function 
obtained from the SDE for constructing the effective charge. Notice the absence of any flavour 
dependence below the 1 GeV region. 

/3-function coefficients, A QCD and (A 2 )), combine in such a way that, when the RG invariant 
combination r(q 2 ) is formed, no net flavour dependence survives in the IR . 

Since, modulo an overall dimensionful factor to render it dimensionless, r(q 2 ) coincides 
with the effective coupling, the origin of this independence can be understood by recalling 
the reason for the Nf dependence of the running coupling in the UV (it should also be noticed 
that the invariant combination r(q 2 ) is related with the UV coupling denned through the 
ghost-gluon vertex in Taylor scheme by nothing but a factor q 2 ). In this case the bigger 
the physical momenta q 2 , the more channels open up for the production of quark anti-quark 
virtual pairs (so that every time a channels opens, the coupling receives a "kick" and goes 
up). However, as soon as q 2 drops below a certain threshold, no energy will be available 
to produce any virtual pairs (not even gluons when q^ < 4ml) so that the residual running 
of the coupling below this value is completely dominated by the IR mass scale introduced 
when defining the effective charge. 

Coming to this specific point, it turns out that [21] one can construct from r(q 2 ) the di- 
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FIG. 7: The effective charge a(q 2 ) defined in Eq. (5.3) for two different values of the IR gluon 
mass: niQ = 500 MeV (solid symbols) and mo = 600 MeV (open symbols). The gray band in the 
background (meant to guide the eye) has been obtained from a continuum extrapolation of the 
N f = 2 + 1 + 1 data. 

mensionless effective coupling a(q 2 ) by pulling out the inverse propagator factor q 2 + m 2 (q 2 ), 
i.e., 

a(q 2 ) = [q 2 + m 2 (q 2 )]r(q 2 ), (5.3) 

which leads to an IR saturating coupling. Notice that this definition is valid for both 
massive and the (already ruled out) scaling solutions (in which case one would have to set 
m 2 (q 2 ) = 0); since in the latter case A(0) — > and F(0) — > oo, the effective coupling does 
not distinguish between the two solutions 6 . 

As a last step, we need to specify the q 2 running of the dynamical mass m 2 (q 2 ). We 
will consider here the simplified setting of [2H E3] under which the mass obeys a power law 



6 As explained in detail in [24 , in the presence of an IR saturating propagator, one should not insist in 
pulling out in front of the effective coupling a simple q 2 factor. Otherwise, one would end up with a 
completely unphysical coupling, namely the one that vanishes in the IR, where QCD is supposed to be a 
strongly coupled theory. 
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running 



m 



m (l) = Tt4; ™o = m(0), (5.4) 

-yZ I ivy) ^ 



and for mo one considers the representative values m = 500 — 600 MeV, consistent with a 
variety of phenomenological studies. The resulting effective charge is plotted in Fig. [7} 

We hasten to emphasize that this is only a toy model, and one should take into account 
that in reality mo differs for different Nf, as clearly seen in the top panel of Fig. [TJ However, 
inserting directly in Eq. (5.4) the saturation value m = A~ 1 (0,/i 2 ) obtained from our sim- 
ulations for different number of flavours Nf, breaks the RG invariance in zero of Eq. (5.3 ) 7 . 
We are evidently in need of better tools for extracting reliable (RG invariant) information 
about the saturation point of the coupling (and probably more data as well in the low 
momentum region); this issue clearly deserves a separate study. 



VI. CONCLUSIONS 



In this paper, we have carried out a systematic and comprehensive analysis of the gluon 
and ghost two-point functions in (Landau gauge) full lattice QCD. 

The configurations used include two light and two light plus two heavy twisted mass 
fermions with masses between 20—50 [MeV] for the light quarks, 95 [MeV] for the strange 
quark and 1.51 [GeV] for the charm quark (in MS scheme at a renormalization scale of 2 
[GeV]). The mass of the lightest pseudoscalar turns out to be between the range of 270 and 
510 [MeV]. As this value does not lie too far from the physical pion mass, it increases our 
confidence in the flavor physics effects reported in this article. Moreover, simulations on 
lattices with up to 48 3 x 96 points, with (3 = 3.90 and 4.20 for N f = 2 and = 1.90 and 
1.95 for Nf = 2 + 1 + 1, allow us to reach momenta down to q ~ 300 [MeV], keeping the 
volume effects under control. 

Our analysis demonstrates that in the intermediate and low momentum region, the gluon 
propagator lessens with the increase in the number of dynamical quarks, whereas, the ghost 
dressing function is enhanced, albeit only slightly. In addition, the heavier a species of 
fermions, the smaller in extent is its effect on the suppression of the gluon propagator. With 

7 One could in principle choose different phenomenological values of mo for different values of Nf but the 
result would not be that different from what is seen in Fig. [7) 
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a heavier enough mass, which prevents its virtual pair production, the fermion fails to screen 
the interaction and gets decoupled from the gluon dynamics altogether. 

When all the pieces of data are put together to construct the effective QCD running 
coupling, a(q 2 ), we observe the behavior anticipated from the massive decoupling solutions, 
namely, a monotonic approach to an IR fixed point. Furthermore, we find that below q ~ 1 
[GeV], this quantity is not directly affected by the variation in the number of dynamical 
fermion flavours. However, considering that an IR gluon mass is introduced while defining 
the effective running coupling, there is indirect dependence on Nf via this mass scale. 

Making the most of the lattice results for the gluon and ghost propagators, we present a 
self-consistent analysis of the ghost 2-point function and extract the unquenched ghost-gluon 



form factor H±. This is a genuine prediction of the SDE study presented in Sect. IV which 
should be confirmed (or refuted) through direct lattice studies of the ghost-gluon vertex. 

Though the data presented here have been obtained for an arbitrary gauge copy selected 
through a gauge fixing algorithm using a combination of over-relaxation and Fourier accel- 
eration, we do not expect that the presence of Gribov copies to alter the conclusions in any 
significant way, given also that their effect is expected to weaken significantly for physical 
volumes as large as the ones considered here. 
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